Properties of trapped neutrons interacting with realistic nuclear Hamiltonians 
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We calculate properties of neutron drops in external potentials using both quantum Monte Carlo 
and no-core full configuration techniques. The properties of the external wells are varied to examine 
different density profiles. We compare neutron drop results given by a selection of nuclear Hamil- 
tonians. including realistic two-body interactions as well as several three-body forces. We compute 
a range of properties for the neutron drops: ground-state energies, spin-orbit splittings, excitation 
energies, radial densities and rms radii. We compare the equations of state for neutron matter for 
several of these Hamiltonians. Our results can be used as benchmarks to test other many-body 
techniques, and to constrain properties of energy-density functionals. 

PACS numbers: 21.30.-x, 21.60.-n, 21.60.De, 21.65.Cd 



I. INTRODUCTION 

There are three major motivations for investigating 
pure neutron systems with external fields ("neutron 
drops") using ab initio approaches. First, neutron drops 
provide a very simple model of neutron-rich nuclei in 
which the core is approximated as an external well acting 
on valence neutrons pJ-0]]. Second, ab initio solutions for 
neutrons trapped by an external potential can be used as 
data for calibrating model effective Hamiltonians and En- 
ergy Density Functionals (EDFs) [jj @ , particularly for 
very neutron-rich systems as occur in astrophysical en- 
vironments like neutron stars. Third, these results may 
serve as useful benchmarks for testing other many-body 
methods. 

These motivations are further supported by the advent 
of new experimental facilities to probe the extremes of 
neutron-rich nuclei, to map out the neutron drip line and 
to inform models of nuclear astrophysical processes 0j- 
Traditional energy density functionals Q are obtained 
by fitting measured properties of stable and near-stable 
nuclei. Extrapolations using these traditional EDFs to 
regions of extreme isospin are sensitive to their less con- 
trolled features that result in large variations in the pre- 
dictions. Beyond these experimental vistas, we desire 
control over properties of EDFs for low-density neutron 
systems, since those properties are important for the pro- 
cesses in the inner crust of neutron stars. It is our long- 
term aim to provide ab initio calculations of trapped neu- 
trons that address these motivations. We report here the 
results with currently available approaches that serve as 
a basis for long-term efforts that will employ improved 
microscopic Hamiltonians and many-body methods. 

We adopt two nucleon-nucleon (AW) interactions 
which fit scattering data and deuteron properties with 
high accuracy, namely the local Argonne v' s (AV8') @ 



and the nonlocal JISP16 [lfj. As shown by accurate cal- 
culations 0, [nHlH , local WW interactions are not suffi- 
cient to describe accurately the properties of light nuclei. 
Even the ground-state of the simplest three-body prob- 
lem, the triton, is significantly undcrbound. 

Different models of three-nucleon interactions (TNIs) 
have been proposed to build a non-relativistic Hamil- 
tonian that reproduces experimental results such as 
ground-state energies, density profiles, and rms radii of 
light nuclei [l4j. TNIs from meson-exchange theory are 
modeled through an operator structure with parameters 
that are fit to experimental nuclear energies. The Ur- 
bana IX TNI (UIX) was fit to 3 H and nuclear matter 
saturation, but it typically underbinds heavier nuclei . 
Other TNI forms, namely Illinois forces, are fit to light 
nuclei 14 1. The most recent is the Illinois-7 (IL7) [19| 
which reproduces nuclear energies up to A = 12 with an 
rms error of 600 kcV. However, the three-pion rings in- 
cluded in IL7 give a strong overbinding of pure neutron- 
matter [2(|. In light of the different data selected for 
tuning these TNIs it is useful to observe their similarities 
and differences with neutron drops. 

JISP16 is a phenomcnological nonlocal WW-intcraction 
written as a finite matrix in a harmonic oscillator (HO) 
basis for each of the WW partial waves. It is constructed 
to reproduce the available WW scattering data using the 
J-matrix inverse scattering approach. In addition, phase- 
equivalent transformations have been used to modify its 
off-shell properties in order to achieve a good description 
of selected states in light nuclei [l(|. It gives a good 
description of most narrow states in light nuclei up to 
about A = 12 [2l| - |23j . However it tends to overbind 
heavier N = Z nuclei ( 16 is overbound by about 15%), 
but tends to underbind as one moves away from the valley 
of stability. 

In this paper we analyze the ground-states and several 
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excited states of neutron drops for four Hamiltonians: 
AV8' without TNI, AV8'+UIX, AV8'+IL7, and JISP16. 
We examine possible sub-shell closures and the spin-orbit 
splittings of odd systems near closed HO shells. We also 
compare results for the neutron matter equations of state 
using AV8' with and without TNIs that could be useful to 
calibrate bulk and gradient terms of Skyrme forces. The 
techniques we use are based on two quantum Monte Carlo 
(QMC) techniques for the Argonne interactions and on 
the No-Core Full Configuration (NCFC) for the nonlocal 
interaction JISP16. We provide quantified uncertainties 
where feasible. 



The Green's function Monte Carlo (GFMC) provides 
accurate energies, radii and and other properties of nu- 
clei up to A — 12 with the Argonne interactions fl9| : 
currently it can be used for systems of up to 16 neu- 
trons. The Auxiliary Field Diffusion Monte Carlo 
(AFDMC) [H[ has similar statistical accuracy as GFMC 
in computing energies of systems of neutrons and can be 
implemented for larger systems, up to more than 100 neu- 
trons [1,[25|]. A comparison between AFDMC and GFMC 
results (obtained with the same Hamiltonian) suggests 
that the systematic uncertainties are of the order of a 
few percent. Improving the trial wave function used in 
AFDMC to include pairing, for example, could further 
reduce these differences. 



For JISP16 we expand the neutron drop wave func- 
tions in a HO basis. For any finite truncation of the 
basis, this provides us with a strict upper limit on the 
total energy of the system. Exact results are obtained by 
considering the limit of a complete (infinite-dimensional) 
HO basis — which we refer to as No-Core Full Config- 
uration (NCFC) [2l|. We can obtain the total energies 
for systems up to A = 14 nucleons to within a percent 
by a simple extrapolation to the complete basis from a 
series of successive finite truncations. The extrapolation 
of other observables such as radii and densities is not 
as straightforward, but for small enough systems we can 
simply consider a large enough basis space in order to 
obtain converged results. Note that in a single run, at 
a fixed truncation, we not only obtain the ground state 
energy, but also the low-lying spectrum of the system. 



II. HAMILTONIANS 

We adopt non-relativistic Hamiltonians with the fol- 
lowing general form: 



-e^+e^^+e^. 

i i<j 



E v *» 

i<j<k 



(1) 

Systems consisting of only neutrons are not expected to 
be self-bound. Therefore it is necessary to include an ex- 
ternal well U ex t(r) in the Hamiltonian to have a confined 
system. We consider both harmonic oscillator (HO) wells 
and a Woods-Saxon (WS) potential. 
The HO wells have the form 



U HO (r) = Kntfr 2 



(2) 



This potential is useful due to its simplicity and the fact 
that the ground state may be driven to arbitrary low 
density (i.e. with arbitrarily weak external harmonic 
potential strength) or arbitrarily large particle number. 
The convergence of both the Monte Carlo and the con- 
figuration interaction methods are improved due to the 
lack of any low-lying states with long-range tails in the 
wave function. This feature enables applying our results 
for tests of EDFs over a range from moderately low to 
rather high densities. Most of our results are for HO 
wells. Woods-Saxon wells have been used in other cal- 
culations of properties of neutron drops. In particular 
neutron drops with a WS well and the Argonne vi& NN 
and Illinois-2 potentials have been shown to provide a 
good description of oxygen isotopes @ . The WS form is 



Uwsir) 



1 _|_ e (r-R)/a ' 



(3) 



where we have used a = 1.1 fm, Uq — —35.5 MeV and 
R = 3 fm, that is, the same parameters as in Ref. [l|. 

In addition to the total energy, E = (H) , we also calcu- 
late other observables, such as the external energy (f7 ext ), 
the internal energy E lnt = (H) — (U ext ), and the rms ra- 
dius r. Note that for a HO well, the external energy 
(f^ext) is proportional to r 2 , and thus the quantities E, 
Emt , and r are not independent observables in a HO well; 
however, in a WS well they are independent. 



A. Argonne NN interaction and three-body forces 



The plan of the paper is the following: in Sec. |TT] wc 
describe various Hamiltonians we consider in this work. 
Then, in Sec. IIIII we briefly review the different Monte 
Carlo many-body techniques used to solve for the neu- 
tron drops. Sec. IIVI presents an overview of the NCFC 
approach and Sec. [V] presents our main results for finite 
neutron drops. We present our results for neutron matter 
in Sec. IVII and our conclusions in Sec. IVIII 



One of the NN potentials we adopt here is the Argonne 
AV8' Slg]. It is a simplified version of the Argonne 
AV18 [23], with the advantage that it can be exactly in- 
cluded in both GFMC and AFDMC algorithms without 
treating any part perturbatively. Other non-local opera- 
tors appearing in AV18 must be included as a perturba- 
tion of AV8' in QMC calculations @ . These perturbative 
corrections can be accurately computed within GFMC, 
but not in AFDMC. Thus we consider AV8' to facilitate 
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comparisons of the two different QMC methods. The 
difference between the binding energies from AV8' and 
AV18 is very small in light nuclei, and much smaller in 
pure neutron systems; about 0.06 MeV per neutron 14 1. 
The Argonne AV8' is a sum of eight operators: 



(4) 



in which the v p {rij) depend on the distance between the 
nucleons i and j, and Ofj are operators. Their form is 



O 



p=i, 



(1, Bi ■ Bj,Sij,Lij ■ Sij) x (1, 



(5) 



with Sij the tensor operator, Lij the relative angular 
momentum and the total spin. The v p (r) parts are 
fit to reproduce the S and P partial waves as well as 
the 3 £>i wave and its coupling to 3 Si of the full AV18 
potential [27 1. 

In this paper we consider the AV8' alone and com- 
bined with two different three-body forces: the Urbana- 
IX (UIX) (HI and the Elinois-7 (IL7) 03- Just like for 
the AW interaction, the TNIs are sums of several opera- 
tors: 



Vi 



— aPW rfiK'PW 

- A 2tt U ijk 



A 2tv U ijk 



A 3tt U ijk 



(6) 



4- A n R -4- /|T=3/2 n fl,T=3/2 
+ A RPi]k + A R U ijk 

Both TNIs include the Fujita-Miyazawa operator (3 27r < pw 
and a phenomenological part O r , while only IL7 has the 
°ijk Pijk > and °i/k terms. In the Fujita- 
Miyazawa term [29[ two pions are exchanged between 
the three nucleons with the creation of an intermediate 
excited state. The phenomenological part O r has no spin 
or isospin dependence. The additional IL7 terms involve 
exchanges of two or three pions and a pure T = 3/2 
repulsion. A full description of the operators is given in 
Refs. 03, Gl]. 

The A™ and An parameters of UIX are determined 
by reproducing the binding energy of 3 H and nuclear 
matter (28[. The UIX model has been used to inves- 
tigate properties of neutron matter (see for example 
Refs. [25|, [3(J [H| and references therein). The resulting 
equation of state will support neutron stars larger than 
two solar masses. 

The Illinois forces are more sophisticated than UIX. 
The A coefficients are determined by fits to binding en- 
ergies of light nuclei 0, EH . The Illinois forces give a 
good description of properties of nuclei up to A = 12, 
including both ground states and excited states, however 
three-pion operators are very attractive in pure neutron 
systems and they overbind neutron matter at large den- 
sities [2(|. In this work, we consider IL7 as described 
in El. 



B. JISP16 

The JISP16 AW interaction is determined by inverse 
scattering techniques from the np phase shifts and is, 



therefore, charge symmetric. JISP16 is available in a 
relative HO basis pi| and can be written as a sum over 
partial waves 



V 



Vs -j- t 12 i nZ > A 



S,J,T 
nl,n'i> 



[n't'l 



(7) 



S,J,T 



where Ml = 40 MeV and J = £ + s. The HO relative 
coordinate wave function is written (r\n£) — lZ n i(r). A 
small number of coefficients {A^ n T f ,} are sufficient to 
describe the phase shifts in each partial wave. Note that 
the JISP16 interaction is non-local and its off-shell prop- 
erties have been tuned by phase-shift equivalent transfor- 
mations to produce good properties of light nuclei. For 
example, JISP16 is tuned in the 3 5i — 3 Di channel to give 
a high precision description of the deuteron's properties. 
Other channels arc tuned to provide good descriptions of 
3 H binding, the lo w-ly ing spectra of 6 Li and the bind- 
ing energy of 16 [lj|. After its initial introduction, it 
was realized that the 16 O energy was not fully converged 
and JISP16 overbinds 16 by about 15 to 20 MeV |l|. 
With these off-shell tunings to nuclei with A > 3 one may 
view JISP16 as simulating, to some approximation, what 
would appear as NNN interaction contributions (as well 
as higher-body interactions) in alternative formulations 
of the nuclear Hamiltonian. 



III. QUANTUM MONTE CARLO METHODS 

Both of our QMC methods use diffusion Monte Carlo 
to project the lowest-energy eigenstate out of a trial wave 
function by a propagation in imaginary time r: 



*(r) 



_ -{H-E T )r 



(8) 



where Et is a normalization factor. In the r — > oo limit 
the only component of that survives is the lowest- 
energy one not orthogonal to ^t- 



*o = hm *(r) . 



(9) 



The evolution in imaginary time is performed by using 
Monte Carlo integration to evaluate 



$(R,t) = J dR'G(R,R',T)^ T (R') . 



(10) 



where G(R, R' , r) is an approximation of the many-body 
Green's function of the Hamiltonian, and R and R' are 
the positions of all A nucleons: R = (ri, . . . , fjv). The 
exact form of G(R, R' , r) is unknown, but it can be accu- 
rately approximated as a product of many G(R, R', At) 
for a small time step, At. The main difference between 
GFMC and AFDMC is in their representations of * and 
the structure of the initial ^t- 
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A. GFMC method and trial wavefunction 

GFMC uses a complete spin-isospin representation of 
the many-body wave function; ^(R,t) is written as a 
vector with 2 A N(T) complex components. Here the 2 A 
allows for all possible nucleon spin up or down combina- 
tions and N(T) is the number of proton- neutron combi- 
nations with the desired isospin. In the case of neutron 
drops N(T) = 1. The exponential growth of the vector 
size with the number of neutrons currently limits GFMC 
calulations of neutron systems to N < 16. 

Calculations of nuclei with realistic interactions face a 
sign problem eventually as the trial wave function is prop- 
agated to large r. To deal with this problem, we use the 
constrained path algorithm to obtain configurations with 
the largest possible overlap with the ground state [32| . 
This method is similar to the fixed node approximation 
in that it is exact in the limit of an exact constraint, and 
it is stable to large imaginary time; however it does not 
provide an upper bound. We then extend the propaga- 
tion without constraint for as long as possible to obtain 
the energy and ground-state properties. The constraint 
and the convergence properties of GFMC are discussed 
in [32| . Figure 3 in that paper shows some convergence 
results for neutron drops. 

The trial wave functions used in our GFMC calcula- 
tions for neutron drops in HO wells arc somewhat sim- 
plified from the ones described in Ref. @ for nuclei: 



|*t) = 



\*j) 



%<3 



\*J) 



\$ N (JMTT 3 )) 



(11) 



(12) 



The non-central Uij and associated central f c are optimal 
correlations for neutron matter of the form described in 
Ref. [33|. For drops with N < 8, the $jv is expanded 
in an LS basis of s- and p-shell oscillator functions as 
described in Ref. @. For N > 8 we use a "BCS" ansatz 
$bcs of the form introduced in Refs. [11 Hi]. For N = 8 
the two forms give very similar variational energies. 

The BCS pairing is important, particularly for low- 
density systems and when trying to calculate even-odd 
staggering of the energies. In this paper we consider 
^bcs for only J = states or, in the case of odd N, for 
states in which the total J is carried by a single neutron. 
Such 4>bcs are written using correlated pairs of neutrons 
with total spin S = and a L = spatial wave function 
4>ij expanded in Os, Op, Is, and Od single-neutron wave 
functions: 

4>ij = Pos<i>Os{ri)<j>os{rj) + Pop<fop(n)<fop(rj)Pi(fij) 
+ Pu<l>u(n)<i>u{rj) 

+ fhd<fod(n)<fod(rj)Pst(rij) ■ (13) 
The /3 n i are variational parameters; only their ratios are 




FIG. 1. (color online) Contours of VMC energies for 12 neu- 
trons in a 10-MeV HO well with AV8'+UIX versus the Od and 
Is coefficients in the BCS wave function. The dots show the 
cases that were computed using GFMC. 



relavant. For even N, 
$bcs(J = 0) 



E 

p 



10*: 



SlS 2 



■ s N ) 



(14) 



where the sum is over all partitions of N neutrons into 
N/2 spin-up neutrons and N/2 spin-down neutrons, i 
is from the set of spin-up neutrons and j is from the 
set of spin-down neutrons, \4>ij\ is a determinant, and 
\s\S2 ■ ■ ■ sn) is the spin-vector component with the given 
spins. For odd N and Sx/2 or d 5 / 2 states we have 

$BCS(M=J) = A[(f>L,S,J,Al{rN,Cr N ) $BCS,N-l(J= 0)] , 

(15) 

which is also a sum over partitions of determinants. 

Figure [T] shows the variational Monte Carlo (VMC) en- 
ergies computed for systems of 12 neutrons with different 
BCS parameters for the Od and Is pairs; the (3q s and /3q p 
are fixed at 100 each so the 8 n core is almost full. The 
optimized choice is used for the GFMC calculations. The 
total VMC or GFMC energy is not very sensitive to the 
choice of parameters (the contour interval is only 0.08%) 
- the effects on pairing can be larger. A detailed descrip- 
tion of the algorithm as well as the importance sampling 
technique (constrained propagation) used to reduce the 
variance can be found in Refs. 0, . 

The VMC energies computed with or even just 4" j 
are much closer to the final GFMC energies than is the 
case for real nuclei. The values using tyj are typically 
less than 10% above the GFMC values and often display 
better convergence as the number of unconstrained steps 
(see Ref. [32|) is increased. 



B. GFMC numerical convergence and error 
estimate 

QMC calculations have an easily quantified statistical 
error arising from the Monte Carlo method. It is not so 
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straightforward to determine the magnitude of the sys- 
tematic errors arising from approximations made in the 
constrained-path implementation. These have been ex- 
tensively discussed for GFMC in Refs. [|,|32|. 

Quantities other than the energy are usually evaluated 
by combining mixed and variational estimates: 



in giving the spinor components, namely 



(* |0|*o) = 2(* |C|*t) - (*t|C|*t}- 



(16) 



No extrapolation for the energy is required since the 
ground state is an eigenstate of the Hamiltonian and the 
propagation commutes with the Hamiltonian. This linear 
extrapolation is usually very accurate if the calculation 
begins with a good trial wave function. It is also pos- 
sible to estimate other obscrvables by forward walking 
techniques or by adding a small perturbation eO to the 
Hamiltonian and computing the difference between the 
energy of the original and perturbed Hamiltonian. We 
have not pursued these methods for the present paper. 
Hence the internal energy results as well as other quanti- 
ties reported below are extrapolated and potentially not 
as accurate as the full energy. 



C. AFDMC method and trial wave function 

The presence of spin operators in the Hamiltonian re- 
quires a summation of all possible good spin states in 
the wave function. In AFDMC the spin states are sam- 
pled using Monte Carlo techniques [24j . This sampling is 
performed by reducing the quadratic dependence of spin 
operators in the exponential of Eq. (|5J| to a linear form by 
means of a Hubbard-Stratonovich transformation. The 
effect of an exponential of a linear combination of spin 
operators consists of a rotation of the spinor for each neu- 
tron during the propagation, and this permits the use of 
a much simpler basis in the trial wave function. The re- 
sult is that the wave function is less accurate, but it can 
be rapidly evaluated. Since both positions and spins can 
be sampled, the AFDMC method can be used to solve for 
the ground state of much larger systems than GFMC - 
more than one-hundred neutrons currently may be solved 
with AFDMC. 

More detailed explanations of the AFDMC method 
and how to include the full NN interactions and TNIs in 
the propagator can be found in Refs. [2(| HH, HU, where 
the constrained-path approximation used to control the 
fermion sign problem is also discussed. 

The AFDMC method projects out the lowest energy 
state with the same symmetry as the trial wave func- 
tion from which the projection is started. The trial wave 
function used in the AFDMC algorithm has the form: 



\* T (R,S)) 



n 



K 3 



\$ JM (R,S)}, (17) 



= «i|t> + dil4.> , 



(18) 



where Ui and di are complex numbers. The Jastrow func- 
tion f(r) is the solution of a Schrodingcr-likc equation for 
f(r < d), 

-—V 2 f(r)+ a v c (r)f(r)=\f(r), (19) 
m 

where v c (r) is the spin- independent part of the nucleon- 
nucleon interaction, and the healing distance d and a are 
variational parameters. For distances r > d we impose 
f(r) = 1. The Jastrow part of the function in our case 
has the only role of reducing the overlap of nucleons, 
therefore reducing the energy variance. 

The antisymmetric part of the wave function is 



$ JM (R,S) = [%2D{<l> a (fi,8i)} 



J.M 



(20) 



where a = {n, I, j, rrij} is the set of quantum numbers of 
single-particle orbitals, and the summation of determi- 
nants gives a trial wave function that is an eigenstate of 
J 2 and M. The single-particle basis is given by 



(21) 



The radial components <& n ij & re obtained by solving the 
Hartree-Fock problem with the Skyrme force SKM [37| . 
Y^ mi are spherical harmonics, and £, s ,m s are spinors in 
the usual up-down basis. For each (J, M) set of quan- 
tum numbers there are several combinations of single- 
particle orbitals. We typically perform several simula- 
tions to identify the ground-state and order of excited 
states. It is also possible to include a BCS pairing term 
in the trial wave function in AFDMC, as has been done 
for GFMC. This would allow more accurate treatment 
of pairing in very low-density systems, and is currently 
under development. 

For neutron matter calculations we change the anti- 
symmetric part of the wave function to be the ground 
state of the Fermi gas, built from a set of plane waves. 
The infinite uniform system is simulated by a cubic pe- 
riodic box of volume L 3 according to the density of the 
system. The momentum vectors in this box are 



(22) 



where a labels the quantum state and n x , n y and n z are 
integer numbers describing the state. The single-particle 
orbitals are given by 



(23) 



where S = (si, . . . , sat). The spin assignments Si consist 



Again, it is possible to generalize the neutron matter cal- 
culations to include BCS pairing in the trial state [38l[39j . 
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D. AFDMC numerical convergence and error 
estimate 

AFDMC is very similar in concept to GFMC and 
convergence and error estimates are also similar. The 
statistical errors are easily evaluated and controllable. 
AFDMC also uses a constrained-path method to circum- 
vent the fermion sign problem, and hence the results de- 
pend to some degree on the choice of trial function. 

Although it is possible to do some unconstrained prop- 
agation with AFDMC, it is more limited than GFMC be- 
cause the sampling of the spins introduces a fermion sign 
problem earlier. Uncertainties can also be addressed by 
choosing several different initial trial states. We hnd the 
constrained path results to be reasonably accurate when 
compared with GFMC for small systems. 

In Table U we compare the GFMC and AFDMC total 
energy results for neutron drops in 5 and 10 MeV HO 
wells. Overall the agreement is of the order of a few 
percent or better. Statistical errors due to sampling can 
be made arbitrarily small, these are about 0.2% or less 
for the AFDMC and GFMC calculations. 

For the lower densities generated by the 5 MeV well the 
GFMC energies are significantly lower than the AFDMC 
energies, by up to 3%. A plausible explanation for these 
discrepancies is the fact that we have not incorporated 
BCS pairing in the AFDMC calculation. Indeed, the dif- 
ferences are nearly zero for the closed shell at N = 8 
(where pairing does not play a role) and grow as we go 
towards the middle of the shell, N = 14. We are there- 
fore pursuing the inclusion of BCS pairing in the AFDMC 
calculations of neutron drops in order to improve the re- 
sults at low densities. On the other hand, for the 10 MeV 
well the GFMC and AFDMC results are all within 1% 
of each other, with the AFDMC typically lower than the 
GFMC energies. 

Systematic errors in calculations of neutron matter are 
similar in spirit. The trial wave function can affect results 
for the energy at low densities where pairing is important. 
At larger densities, though, pairing provides a very small 
fraction of the total energy of the system and calculations 
are much less sensitive to the choice of the trial state. 

In addition to the Monte Carlo errors and the depen- 
dence on the choice of the trial state, we have to con- 
sider finite-size effects for calculations of neutron mat- 
ter. We enforce periodic boundary conditions and fix 
the number of neutrons to be a closed shell in the pe- 
riodic free-particle basis. In order to reduce finite-size 
effects we performed simulations with 66 neutrons. Free 
fcrmions for TV = 66 provide a kinetic energy very close 
to the infinite limit. Any possible finite-size effect due to 
the truncation of the potential energy is properly taken 
into account by considering several replicas of the simu- 
lation box as described in Ref. [2(| HH . Typically these 
uncertainties are very small for bulk properties like the 
ground-state energy. 
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14 


0+ 


142.2(2) 
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TV 


,r 


10 MeV HO well 



3 
3 
4 
5 
5 
6 
7 
7 
8 
9 
9 
10 
11 
11 
12 
13 
13 
14 



TABLE I. Comparison of GFMC and AFDMC total ener- 
gies for neutron drops in 5 MeV and 10 MeV wells with 
AV8'+UIX. Statistical errors due to the Monte Carlo sam- 
pling are given in brackets. 



IV. NO CORE FULL CONFIGURATION 
METHOD 

The NCFC method is based on a series of no-core con- 
figuration interaction calculations with increasing basis 
dimensions. In this approach the wave function for the 
V neutrons is expanded in an TV-body basis of Slater 
determinants of single-particle states, and the many- 
body Schrodingcr equation becomes a large sparse matrix 
eigenvalue problem. We obtain the lowest eigenstates of 
this matrix iteratively. In a complete basis, this method 
would give exact results for a given input interaction V. 
However, practical calculations can only be done in a 
finite-dimensional truncation of a complete basis. We 
perform a series of calculations until we reach numeri- 
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cal convergence in a sufficiently large basis space, or wc 
employ a simple extrapolation [2l| to the complete basis. 



A. Description of basis space 

Our choice for the basis is the harmonic oscillator (HO) 
basis so there are two basis space parameters, the HO en- 
ergy huj and the many-body basis space cutoff iV max . The 
cutoff parameter N max is defined as the maximum num- 
ber of total oscillator quanta allowed in the many-body 
basis space above the minimum for that number of neu- 
trons. Numerical convergence is defined as independence 
of both basis space parameters iV max and fku, within eval- 
uated uncertainties. Note that the basis space parameter 
fku is not necessarily the same as that of the HO well Ml 
that confines the neutrons. 

We employ a many-body basis in the so-called M- 
scheme: the many-body basis states are Slater determi- 
nants in a HO basis, limited by the imposed symmetries 
— parity and total angular momentum projection (M), 
as well as by N max . Each single-particle HO state has its 
orbital and spin angular momenta coupled to good total 
angular momentum, j, and magnetic projection, m. Here 
we only consider natural-parity states, and utilize M = 
for an even number of neutrons, and M = \ for an odd 
number of neutrons. In this scheme a single calculation 
gives the entire spectrum for that parity and N max . 

The NCFC approach satisfies the variational principle 
and guarantees uniform convergence from above to the 
exact eigenvalue with increasing iV max . That is, the re- 
sults for the energy of the lowest state of each spin and 
parity, at any iV max truncation, are strict upper bounds 
on the exact converged answers and the convergence is 
monotonic with increasing iV max . 

The challenge for this approach is that the matrix 
dimension grows nearly exponentially with increasing 
N max . The calculations presented here have been per- 
formed with the code MFDn j40T - l42l ] which has been 
demonstrated to scale to over 200,000 cores. For small 
neutron drops we are able to achieve converged results to 
within a fraction of a percent by using a sufficiently large 
basis space, at least for neutrons in a HO well of 10 MeV 
and above. In order to achieve converged NCFC results 
directly (i.e. without extrapolation) for more than 10 
neutrons using JISP16, we would need to obtain eigen- 
states of matrices that are beyond the reach of present 
technologies. However, for up to 22 neutrons, we can 
utilize a sequence of results obtained with N max values 
that are currently accessible, in order to extrapolate to 
the infinite or complete basis space limit. 



B. NCFC numerical convergence and error 
estimate 

We carefully investigate the dependence of the results 
on the basis space parameters, iV max and fko. Our goal is 




_i i i i i i i i i 

10 15 20 25 30 
basis space hco (MeV) 



FIG. 2. (color online) Ground state energy for 10 neutrons 
in a HO well of 10 MeV (top) and 20 MeV (bottom) for a 
series of finite basis space calculations with JISP16. Note that 
our results in the 20 MeV well converge much more rapidly 
than in the 10 MeV well. The symbols represent extrapolated 
results as discussed in the text and the error bars signify our 
uncertainty estimate at that basis space hw; the yellow band 
represents our final result including numerical error estimates. 



to achieve independence of both of these parameters as 
that is a signal for convergence — the result that would 
be obtained from solving the same problem in a complete 
basis. For the total energy, (H), the guarantee of mono- 
tonic convergence from above to the exact total energy 
facilitates our choice of extrapolating function. 

We use an extrapolation method that was found to be 
reliable in light nuclei: a constant plus an exponential in 
^Vniax [HI, Hj|. That is, for each set of three successive 
7V max values at fixed Ml, we fit the ground state energy 
with three adjustable parameters using the relation 

E gs {N max ) = aexp(-cJV max ) + E gs {oo) . (24) 

Under the assumption that the convergence is indeed ex- 
ponential, such an extrapolation should get more accu- 
rate as iV max increases; we use the difference between the 
extrapolated results from two consecutive sets of three 
N max values as an estimate of the numerical uncertainty 
associated with the extrapolation. 

For a reasonable range of basis space parameters Hco, 
this assumption appears to be valid, as is evident from 
Fig. [3 In this figure we show the ground state energies 
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basis space h(0 (MeV) basis space hco (MeV) 



FIG. 3. (color online) The rms radius (left) and internal energy (right) for 10 neutrons in a HO well of 10 MeV (top) and 20 
MeV (bottom) for a series of finite basis space calculations with JISP16. Note that our results in the 20 MeV well converge 
much more rapidly than in the 10 MeV well. The yellow band represents our best estimate (including an error estimate) for 
the infinite basis space result. 



for 10 neutrons in a 10 MeV and 20 MeV HO well for a 
series of finite bases, as well as the extrapolations with 
their uncertainties indicated by error bars. The error bars 
on the extrapolations using calculations up to 7V max = 10 
are all smaller than the corresponding error bars from 
calculations up to 7V max = 8; and the error bars on the 
extrapolations using calculations up to iV max = 12 are all 
smaller than the corresponding error bars from iV max = 
10. Furthermore, we see that the dependence on the basis 
space Huj decreases with increasing -/V max , and that the 
extrapolated results for a given A^ max agree within each 
other's error bars. Our total error estimate is based on 
a 5 MeV region in hu> which has the smallest error bars 
and minimal fko dependence. 

In order to perform this extrapolation to the infinite 
basis space, we need finite basis space calculations up to 
]V max = 8 or higher. Above 22 neutrons, our calculations 
are limited to N max = 4, so we only have variational 
upper bounds for the total energy of these larger neutron 
drops. For the 10 MeV HO well, our results are not yet 
converged at this basis space, but for the 20 MeV HO 
well, these upper bounds are likely to be within a few 
percent of the converged results. 

The internal energies and the rms radii do not converge 
monotonically with A^ max , in contrast to the total energy. 
Currently, we do not have a reliable method to perform 
the extrapolation to the infinite basis space for these ob- 



servables. The rms radius seems to converge from above 
for small basis space parameters Huj, but from below for 
large basis space parameters, see the left panels of Fig. EI 
Hence there is a 'sweet spot' in the basis space parameter 
fko for which the radius and equivalently, the external en- 
ergy (fext), is approximately independent of iV max . We 
use our results in this region as an estimate of the infinite 
basis space result, with error bars based on the residual 
N max and hjj dependence in a window around the 'sweet 
spot', as indicated by the yellow band in the left panels 
of Fig. El 

The internal energy appears to converge from above, 
at least for hu> values in the region that is optimal for 
the total and external energies. However, we have not 
been able to find a robust convergence pattern; e.g. an 
exponential extrapolation does not work very well, as 
can be seen in the top right panel of Fig. [3l The most 
reliable estimate for the internal energy in the infinite 
basis space appears to be the difference between the (ex- 
trapolated) total energy and the external energy based 
on the convergence at the 'sweet spot' explained above, 
E-mt = (H) — (t/ cx t). This is depicted by the yellow band 
in the right panels of Fig. [3J 

In order to make a meaningful estimate of the con- 
verged (NCFC) results for the rms radius and the exter- 
nal and internal energies, we need to perform a set of 
calculations for a range of basis space parameters fku at 
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least up to N max = 10. We therefore give these results 
only up to 14 neutrons. 



V. RESULTS FOR NEUTRON DROPS 
A. Total energy 

In Table Ql] we present the principal results of this 
study: the total energies for neutron drops confined in 5 
MeV, 10 McV, and 20 MeV HO wells with the AV8'+UIX 
and JISP16 potentials. These HO wells are convenient 
for ab-initio calculations because one can probe very low 
to very high densities with a simple asymptotic form of 
the wave function and an arbitrary number of neutrons 
can be bound in the well. In order to provide a very dif- 
ferent probe of density functionals in the extreme isospin 
limit, we also include select results in a WS well with the 
AV8'+UIX potential. 

We show the lowest + energy for even N and low- 
est values for several J™ for odd N . We present results 
only for natural parity states. The AV8'+UIX values 
up to iV=14 were computed by GFMC while the larger 
drops were computed using AFDMC, with the exception 
of the 20 McV HO well results. There are no results 
from GFMC available for the 20 McV HO well, due to 
the strong fcrmion sign problem with that external field 
strength; those results were all obtained using AFDMC. 

The JISP16 values were all computed by NCFC. There 
are no results from NCFC available in the 5 MeV trap 
above 14 neutrons due to poor convergence with available 
computer resources. Above 22 neutrons we only provide 
strict upper bounds; for the 20 MeV HO well we expect 
the converged energies to be within a few percent of these 
upper bounds. 

Figure 0] shows the energies of N neutrons in two dif- 
ferent HO wells, scaled by Ml V< 4 / 3 ); for odd N only the 
lowest energy found is shown. The scaling by Ml A^ 4 / 3 ) 
is motivated by the expected results in local density ap- 
proximation. The factor V 4 / 3 comes from the tradi- 
tional scaling with N times the increase in potential en- 
ergy arising from the increase in radius of the system 
with particle number proportional to V 1 / 3 . In addition 
to the AV8'+UIX and JISP16 values presented in Ta- 
ble HH we also show results for AV8' without any TNI 
and AV8'+IL7. All interactions show a very pronounced 
peak for three neutrons, and dips at the expected HO 
magic numbers, N = 2, 8, 20, and 40. The dips at the 
HO magic numbers are expected due to the HO nature 
of the confining well. 

With an equation of state of the form E = £j—k F 
with k F = [37rV] 1/3 , the energy is given by the Thomas- 
Fermi expression: E TF = ^/ 2 Ml (3V)( 4 / 3 )/4. For free 
fermions (£ = 1) the Thomas-Fermi results would be a 
horizontal line at 3 4 / 3 /4 « 1.081, for a unitary Fermi gas 
with £ = 0.4 the Thomas-Fermi results would be a hori- 
zontal line at 0.684. The calculated results are all below 
the free Fermi gas (even for the case of three neutrons) 
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FIG. 4. (color online) Energy of the lowest neutron drop 
states confined in a HO well with = 10 MeV (top) and 
hQ, = 20 MeV (bottom) as a function of the number of neu- 
trons. Results for AV8' (plus TNI) where obtained using 
AFDMC, with MC statistical error bars as well as a band in- 
dicating the 1% systematic uncertainty discussed in Sec. lIII Dl 
results for JISP16 are obtained from NCFC with error bars 
reflecting the total numerical uncertainty, and strict upper 
bounds obtained with NCSM in finite basis spaces. Note the 
pronounced dips at the expected HO magic numbers N — 2, 
8, and 20. 



since the interaction is attractive. All our results are 
above the unitary Fermi gas because there are significant 
finite-range corrections for neutron matter. In addition 
repulsive gradient terms in the density functional are re- 
quired to reproduce the ab-initio results @. A detailed 
investigation of these effects is being pursued. 

From Fig. g] it is evident that adding UIX to AV8' 
increases the energies of neutron drops, whereas IL7 de- 
creases the energies. These results were expected; the 
two-pion part of UIX is attractive in the isospin T = 1/2 
triples that appear in nuclei. However neutron drops con- 
tain only T = 3/2 triples for which the two-pion part is 
very small [3ll l44j : this leaves only the repulsive cen- 
tral part of UIX. On the other hand, IL7 contains the 
three-pion term that is strongly attractive in T = 3/2 
triplets Q. 

The energies with the nonlocal 2-body interaction 
JISP16 are generally below the AV8' results, but above 
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5 MeV HO well 


10 MeV HO well 


20 MeV HO well 


WS well 


N 


,r 


AV8'+UIX 


JISP16 


AV8'+UIX 


JISP16 


AV8'+UIX 


JISP16 


AV8'4-UIX 


3 


2 


22.89 


22.73(1) 


46.69(1) 


46.512 


97.1(1) 


98.094 




3 


2 


22.61 


22.40(1) 


45.48(0) 


44.833 


91.7(1) 


90.915 




4 


0+ 


29.99 


29.69(1) 


62.04(1) 


60.842 


131.1(1) 


126.31 




5 


2 


41.22(1) 


40.65(15) 


84.02(2) 


82.86(2) 


175.2(1) 


173.00 




5 


2 


41.02 


40.4(2) 


82.97(1) 


80.68(2) 


169.5(1) 


162.71 




6 


0+ 


48.52(1) 


47.6(2) 


98.95(2) 


95.74(3) 


205.8(2) 


193.64 


-80.6 


7 


2 


59.17(1) 


57.9(2) 


118.9(0) 


115.67(5) 


246.4(2) 


237.11 


-90.9(1) 


7 


3 
2 


59.73(1) 


58.5(2) 


121.1(0) 


118.9(1) 


254.7(2) 


249.85 


-88.6(1) 


8 


+ 


67.01(1) 


65.4(3) 


135.8(0) 


132.5(1) 


287.4(2) 


278.32(1) 


-103.9(1) 


9 


1 + 
2 


80.92(4) 


78.9(1.5) 


163.7(1) 


159.6(4) 


349.8(2) 


334.32(1) 


-107.8(1) 


9 


3 + 
2 




80.0(1.5) 




162.8(6) 


354.5(2) 


344.42(1) 




9 


5 + 
2 


81.20(3) 


79.3(1.5) 


163.2(1) 


159.4(4) 


343.9(2) 


331.15(1) 


-106.6(1) 


10 


0+ 


92.14(8) 


90.(1.5) 


188.1(6) 


182.1(5) 


400.5(2) 


380.41(1) 


-113.4(1) 


11 


]_ + 
2 


105.9(1) 




216.1(3) 


208.4(1.0) 




434.38(5) 


-115.9(2) 


11 


3 + 
2 








208.0(1.0) 




430.10(4) 




11 


2 


106.3(1) 




217.0(3) 


207.9(1.0) 




430.41(4) 


-116.9(2) 


12 


0+ 


118.1(1) 


116.(6) 


242.(1) 


230.0(1.0) 


509.1(4) 


477.05(5) 


-123 6(3") 


13 


i -4- 
2 


130.8(1) 




268.0(1.0) 


255.8(1.0) 




529.07(6) 


-125.0(3) 


13 


3 + 
2 








256.3(1.0) 




528.74(6) 




13 


5 + 
2 


131.5(1) 




267.6(6) 


255.7(1.0) 




524.77(6) 


-125 9(3) 


14 
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142.2(2) 


140.(10) 


291.9(2) 


277.5(1.4) 




569.3(1) 


-131.6(7) 


15 


2_ -f 
2 


160.1(1) 




316.3(2) 


303.(5) 




619.4(6) 




15 


3 + 
2 


159.1(2) 




320.2(2) 










15 


5 + 
2 


160.0(1) 




317.0(2) 








-139.3(3) 


16 


0+ 


171.6(1) 




341.5(2) 


326.(6) 


730.3(3) 


667.7(6) 


-142.4(7) 


17 


2_ + 
2 


185.5(2) 




368.8(3) 






725.1(8) 




17 


3 + 
2 


183.9(2) 




366.5(3) 


352.(5) 






-148.8(2) 


17 


5 + 
2 


184.9(2) 




371.1(2) 










18 


+ 


195.6(2) 




392.6(3) 


377.(7) 




781.(1) 


-155.1(4) 


19 


I 4- 
2 


209.4(2) 




420.1(2) 


407.(7) 


919.0(3) 


850.(2) 




19 


2 


208.4(2) 




417.9(3) 


403.(7) 


914.9(4) 


838.(1) 


-159.6(3) 


19 


r 4- 

2 


210.0(3) 




422.1(2) 


408.(8) 


926.7(4) 


855.(2) 




20 


+ 


219.9(3) 




441.7(4) 


430.(10) 


976.0(4) 


894.(2) 


-165.0(1) 


21 


7 
.: 






476.8(4) 


465.(25) 




956.(4) 




22 


0+ 


254.(1) 




510.5(5) 


495.(25) 


1123.3(7) 


1018.(5) 




24 


0+ 


289.1(7) 




578.9(5) 


< 596. 


1268.(1) 


< 1144. 




26 


0+ 


324.0(8) 




645.0(6) 


< 660. 




< 1266. 




28 


0+ 


355.(1) 




707.6(7) 


< 723. 


1551.(1) 


< 1379. 




30 


0+ 


390.0(8) 




776.0(9) 


< 786. 




< 1499. 




32 


0+ 


422.(1) 




843.5(9) 


< 847. 




< 1614. 




34 


0+ 


453.(1) 




909 9(9) 


< 914. 




< 1750. 




36 


0+ 


486.(1) 




982.7(8) 


< 986. 




< 1895. 




38 


+ 


514.(1) 




1046.4(8) 


< 1057. 




< 2037. 




40 


0+ 


546.(1) 




1114.3(9) 


< 1128. 




< 2177. 




42 


0+ 


591.(1) 




1197.1(8) 


< 1206. 




< 2320. 




44 


0+ 






1278.(1) 






< 2473. 





TABLE II. Total energy with AV8'+UIX and JISP16 for several different confining wells. Error bars for the AV8'+UIX results 
are statistical only; in addition we expect systematic errors of a few percent, as discussed in Sec. IIIIDI Error bars for the 
JISP16 results are total error estimates. (Errors that are not shown are less than 1 in the last digit.) 
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the AV8'+IL7 ground state energies. In fact, in the 10 
MeV HO trap, the JISP16 results are nearly identical to 
those with AV8'+IL7 up to about 12 neutrons; as the 
number of neutrons increases, the results with JISP16 
deviate more and more from the AV8'+IL7 results. In 
the 20 MeV HO trap, for which we have more accurate 
results with JISP16, the results with JISP16 and with 
AV8' without TNI arc quite similar, even in the sd-shcll 
and beyond. The trend of the upper bounds obtained 
with JISP16 follows the trend of the AV8' ground state 
energies through the sd-shell and into the p/-shell, both 
in the 10 MeV and in the 20 MeV trap. 

As discussed in Sec. Ill Al and in Sec. |VI] below, recent 
studies of the neutron star mass-radius relationship (3ll . 
EU suggest that, at least at higher densities, the AV8' 
+ UIX interactions gives a reasonable neutron matter 
equation of state. The requirement of a two-solar mass 
neutron star implies a repulsive three-neutron interaction 
at moderate and high densities. 

On the other hand, AV8'+IL7 gives a much better de- 
scription of the ground state energies, spectra, and other 
observables for light nuclei (up to A = 12) than either 
AV8' or AV8'+UIX. This may be why the results with 
AV8'+IL7 and with JISP16 (which also gives a good de- 
scription of light nuclei) are quite similar below 12 neu- 
trons. However, none of these interactions have been fit 
to any data beyond the p-shell, and it is unclear which of 
these interactions is more realistic for the neutron drops 
in the N = 8 to N = 40 range. At larger densities 
AV8'+IL7 is too attractive as discussed below. 

In the 10 MeV well, the dips in the energies at N = 16 
and N = 32 suggest subshell closure with AV8'+IL7, but 
not with AV8'+UIX, while the results for AV8' show a 
hint of subshell closure at N = 32. The IL7 TNI docs 
provide a larger spin-orbit splitting than the UIX three- 
nucleon interaction. Similarly, the energies with JISP16 
suggest subshell closure at N = 16 and N = 32 in the 20 
MeV well. The JISP16 results in the 10 MeV well are not 
quite accurate enough to draw firm conclusion regarding 
subshell closure; and we have insufficient results for AV8' 
in the 20 MeV well. 

Somewhat surprisingly, there is no indication of sub- 
shell closure at N = 28. In other words, these results 
seem to suggest closure of the combined f? and p3 sub- 
shell at TV = 32, rather than closure of just the fx at 
TV = 28. Note that the closure of the combined da and 

2 

si subshell at N = 16 corresponds to the recently dis- 
covered subshell closure at 24 O 1461. 



B. Energy differences 

In Fig. [5] we show the difference in total energy between 
neutron drops with TV and with TV — 1 neutrons. Wc 
clearly see the effect of the HO shells: jumps at 2, 8, 
and 20 neutrons, at which the next neutron has to go 
to the next HO shell. Without interactions between the 
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FIG. 5. (color online) Single energy differences in a 10 MeV 
(top) and 20 MeV (bottom) HO well. Results of different 
Hamiltonians are compared. AFDMC and GFMC error bars 
are statistical only; NCFC error bars reflect the total nu- 
merical uncertainty. Horizontal line segments indicate energy 
differences expected from pure HO energies. 



neutrons, we would still have this shell structure, but 
within each shell, all single energy differences would be 
equal, as indicated by the solid reference lines in Fig. [5] 
That is, the gross feature of shell structure arises from 
the confining well and is evident in the plot of the single 
differences as a jump in the calculated energy differences 
as one goes from one shell to the next. 

The detailed fluctuations within a shell are entirely 
due to the neutron interactions. The most promi- 
nent feature is the neutron pairing, in particular in 
the p-shell and also in the (beginning) of the sd- 
shell. This effect can be seen more clearly by look- 
ing at the double difference in total energy A (A) = 
(-1) N+1 [E(N) - \ (E(N - 1) + E(N + 1))], see Fig. E 
The phase (— 1) N+1 is included to make the pairing posi- 
tive definite in the standard BCS theory Without inter- 
actions, the double differences would be zero, except at 
the magic numbers 2, 8, and 20. 

Overall, the pairing seems to decrease as N increases, 
except for the closed (sub)shclls. Note that the pairing 
in nuclei also decreases for larger nuclei [47[ . On the 
other hand, the numerical uncertainties increase with A, 
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FIG. 6. (color online) Double energy differences A(N) = 
(-1) N+1 [E(N) - \{E{N - 1) + E(N + 1))] in a 10 MeV 
(top) and 20 MeV (bottom) HO well. Results of different 
Hamiltonians are compared. AFDMC and GFMC error bars 
are statistical only; the NCFC error bars are omitted for the 
10 MeV HO well, because they would cover the entire vertical 
range for 12 neutrons and above, though a significant part of 
the NCFC numerical error is systematic, and cancels between 
neighboring neutron drops; for completeness, we did include 
the total numerical uncertainty for the NCFC results in the 
20 MeV HO well. 



preventing us from obtaining meaningful results for the 
pairing beyond 22 neutrons using the NCFC approach. 
AFDMC calculations of the pairing gaps will be more 
reliable once BCS correlations have been included in the 
trial state and this is being pursued. For all methods 
we expect that the error in neighboring neutron drops is 
correlated, resulting in a reduced error in calculations of 
energy differences and pairing. 

Despite the numerical uncertainties, there are some 
features that are likely to be robust in Fig. [51 As ex- 
pected, the peaks in the double difference A(V) at the 
magic numbers 8 and 20 stand out for all of the inter- 
actions for which we have results, in particular for the 
10 MeV HO well. In addition, our results suggest sub- 
shell closure at N = 16 for AV8' without TNIs, with 
AV8'+IL7, and with JISP16, but not with AV8'+UIX. 
This closed subshell corresponds to the recently discov- 
ered subshell closure at 24 O [46[, in which TNIs play a 



crucial role. 

In addition to the closure at N = 16, we see evidence 
for subshell closure at N = 6 (the pa) in the 20 MeV 
HO well both with AV8'+IL7 and with JISP16, but not 
in the 10 MeV HO well. We do not have sufficient data 
yet to examine the expected closure of the fx at = 28, 
which was not evident in the plots of the total energy (see 
Fig. |4]), nor for a more detailed analysis of the closure at 
N = 32 suggested in Fig. g] 



C. Level splittings 

If we look at the single-particle and single-hole states 
at the beginning and end of the p-shell, see Fig. [71 we 
find that the spin-orbit splitting between the \ and | 
increases with Ml for all interactions. For three neutrons, 
the splitting between these levels is almost the same for 



a -4 



-12 



-16 



20 



10 



3 neutrons 

1 -particle states 



IT 



7 neutrons 
1-hole states 



x x JISP16 

. . AV8' + UIX (AFDMC) 
• AV8' + IL7 (AFDMC) 
. .. AV8' (AFDMC) 



_i i i i i 

5 10 15 

external field h£2 (MeV) 



"■5 



a 

20 



> 



< 



-10 



d3/2 - d5/2 




sl/2-d5/2 



19 neutrons 
1-hole states 



d3/2-d5/2"'J. 



10 15 
external field hn (MeV) 



20 



FIG. 7. (color online) Spin-orbit splitting in the p-shell (top) 
and level splittings in the sd-shell (bottom) and as a function 
of external field strength. Results of different Hamiltonians 
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JISP16 and AV8'+IL7; however, for seven neutrons the 
splitting is significantly enhanced with IL7. On the other 
hand, AV8' without TNI and AV8'+UIX have almost the 
same splitting. 

The systematic increase in level splittings with increas- 
ing Ml can be understood as follows: With increased Ml, 
the radial shape is increasingly constrained by the HO 
potential and the associated gaussian falloff of the ra- 
dial densities in the surface region. This increase in level 
splittings with Ml may then be interpreted as a conse- 
quence of the increased density gradient in the surface 
region. 

In the sd-shell the splitting between the ds and s 1 lev- 
els (solid lines in Fig. [7]) is much smaller than the splitting 
between these two levels and the d,3 level, in particular for 
AV8'+IL7. This confirms the subshell closure at N = 16 
that was evident from the pairing, see Fig. [SI It is also in 
apparent agreement with the observation in known nuclei 
that the subshell closure at 16 neutrons (both ds and s i 
levels filled) is much stronger than the subshell closure 
at 14 neutrons (only the ds level filled). 

Furthermore, notice that the level ordering can change 
as the strength of the HO well increases in the case of 9 
neutrons: in the weakest well of 5 MeV (i.e. at very low 
density), the si is slightly below the ds level, but as Ml 
increases, the ds becomes the lowest level. Interestingly, 
this happens both with JISP16 and with AV8'+IL7, 
whereas with AV8'+UIX and with AV8' (without TNIs) 
the ds and si are basically degenerate for the 5 MeV 
HO well. 

In the pf-shell we find qualitatively similar results 
with JISP16: a large spin-orbit splitting between the 
fz and fs levels and between the p3 and pi levels, 
a smaller splitting between the pi and fs levels, and 
an even smaller splitting between the fi and p3 lev- 
els. All of these level splittings increase significantly 
with the strength of the HO well: at Ml = 5 MeV, the 
splittings are almost negligible, less than an MeV, and 
within the numerical uncertainty. On the other hand, at 
Ml = 20 MeV (the largest value that we have considered) 
the spin-orbit splittings are of the order of ten(s) of MeV. 



D. Internal energies and radii 

In Table Hill we list our results for the internal energy 
^int = (H) — (E/cxt), as w cll as for the rms radii for sys- 
tems up to 14 neutrons in a HO well with JISP16 and 
with AV8'+UIX, as well as in a WS well with AV8'+UIX. 
Note that for neutron drops in a HO well the radius is di- 
rectly related to the external energy, (f/ xt) = \muj 2 (r 2 ) 
for a HO external field. Overall, the internal energy is 
typically slightly less than half of the total energy (see 
Table HI1 for comparison), and (t/ xt) is slightly more than 
half the total energy. This is to be expected, since the 
total energy scales approximately as MIN^/ 3 ^ 1 cx p 2 ^ 3 , 
and for all cases where the equation of state is propor- 
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FIG. 8. (color online) Internal energy for up to 14 neutrons 
in a HO trap with AV8'+UIX and with JISP16. For details 
see Table Hn] 
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FIG. 9. (color online) Radii for the lowest J states up to 14 
neutrons in a HO trap with AV8'+UIX and with JISP16. For 
details see Table [Till 



tional to p 2 ^ 3 , the virial theorem will give equal internal 
energies and one-body potential energies, each one-half 
of the total energy. 

In Fig. [5] we show the internal energy E- mt , scaled by 
Ml N^/V, of the lowest J = and J = § states for 
up to 14 neutrons in a HO well. In the 10 MeV trap 
the JISP16 and the AV8'+UIX results are rather close 
to each other, significantly closer than the total energies 
shown in Fig. [4] Apparently, the larger differences ob- 
served in Fig. |4] arise primarily from differences in their 
(f^ext) energy shifts. Indeed, the corresponding rms radii, 
and thus (U cx t), start to deviate from each other above 
N = 10, see Fig. [9l The two interactions also give quite 
similar internal energy results in the 5 MeV trap as seen 
in Fig. given the rather large error bars of the NCFC 
results, and the corresponding radii arc almost identical, 
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TABLE III. Internal energies (in MeV) and rms radii (in fm) of neutron drops with various external potentials and a selected 
set of Hamiltonians. The results are plotted in Figs. [8] and [9] The particular many-body method used to produce these results 
depends on the external field, the number of neutrons and the Hamiltonian as described in the text. 



at least up to 10 neutrons. 

Tabic Unl shows that both the internal energies and the 
rms radii in the WS well are of the same order as those in 
the 10 MeV HO trap, even though the total energies are 
very different. In fact, the rms radii are nearly identical 
for the three p-shell neutron drops (N = 6, 7, and 8), but 
in the sd-shell there are significant differences between 
the HO and the WS radii. 

We note that the internal energy Fig. [5] displays a sim- 
ilar odd-even effect due to pairing as we observed in the 
total energy in Fig. 2J The radii of the J = and J = \ 
states also show an odd-even effect, but only in the p- 
shell, for three to seven neutrons; there is no a significant 
odd-even effect for these states above TV = 8 in Fig. [9] 
Also note that the radii of states with different J in odd 
neutron systems are slightly different, in particular in the 
p-shell, and more so with JISP16 than with AV8' + UIX, 
as can be seen in Table IIIII 



E. Densities 

We present a sample set of radial density distributions 
computed with JISP16 and with the AV8'+UIX Hamilto- 
nian in Fig.[TUJ The band thickness of the JISP16 results, 
obtained with NCFC, is our best estimate of the total nu- 



merical uncertainty in these densities; the AV8'+UIX re- 
sults were calculated with GFMC, and the error bars cor- 
respond to the statistical errors in the GFMC approach. 
Given the HO nature of the trap, all density distribu- 
tions fall like gaussians at distances sufficiently far from 
the origin. 

The various densities for 8 neutrons (top panel of 
Fig. [Hjl closed p-shell) are quite similar for the two differ- 
ent interactions. The only difference is that the central 
densities, below 1 fm, are about 10% to 20% higher with 
JISP16 than with AV8'+UIX, but the shape is essentially 
the same, and above 2 to 3 fm the densities are practi- 
cally on top of each other. This could be expected from 
the similar rms radii for these cases shown in Table IIIII 
As the HO trap strength increases, the density distribu- 
tion gets compressed, the rms radius decreases, and the 
central density increases, as one would expect. The ra- 
dial shape, a slight dip at r = 0, is typical for the closed 
p-shell, and is qualitatively the same for weak and strong 
HO traps. 

However, the shape of the density profile for 14 neu- 
trons is somewhat different for the two interactions; fur- 
thermore, the shape seems to depend on the strength of 
the HO trap, at least for JISP16. With JISP16 in the 20 
MeV trap (dashed curve in the bottom panel of Fig. [TO]), 
the density has a clear dip at the center, and peaks at a 
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FIG. 10. (color online) Radial density distributions for 8 (top) 
and 14 (bottom) neutrons in different HO traps with JISP16 
and with AV8'+UIX. 



distance of about 1 fm from the center. In fact, the shape 
of this density is rather similar to that of 8 neutrons in a 
HO trap. On the other hand, in the 10 MeV trap there is 
no evidence for such a dip; rather, within the estimated 
numerical accuracy, the density seems to fall off mono- 
tonically from a central value of about 0.11 to 0.12 fm -3 
with JISP16. On the other hand, the densities obtained 
with the AV8'+UIX Hamiltonian for 14 neutrons seem 
to be slightly enhanced in the central region: in the 10 
MeV trap the central density with with the AV8'+UIX 
is about 20% higher than with JISP16. 

This difference between the density profiles of 14 neu- 
trons in a HO trap with JISP16 and AV8'+UIX could 
well be related to the presence (JISP16) and absence 
(AV8'+UIX) of sub-shell closure for 16 neutrons; and 
both are likely to be related to differences in spin-orbit 
splittings. It would be interesting to compare these den- 
sities with those obtained with other realistic potentials, 
and in particular to investigate the effect of different 3- 
body forces on the density profiles as well as on the spin- 
orbit splittings and sub-shell closures. 
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FIG. 11. (color online) Equation of state of neutron matter 
as a function of the density for different Hamiltonians. 
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0.16 
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19.10(2) 


10.62(3) 



TABLE IV. Equation of state of neutron matter as a function 
of the density for various Hamiltonians. 



VI. NEUTRON MATTER 

The equation of state of neutron matter is important 
to properly fix the bulk term of Skyrme-type EDFs. We 
report the AFDMC results for the energy per neutron as 
a function of the density in Table ITVl and display them in 
Fig. [TU For the AFDMC Quantum Monte Carlo calcu- 
lations, the results are for a system of 66 neutrons with 
periodic boundary conditions. The calculation is very 
similar to those of the neutron drops, except the single- 
particle orbitals in the trial wave function are replaced 
by plane waves that respect the periodic boundary condi- 
tion, as described in Sec. MI Cl More details can be found 
in Refs. [2(1 [2f|. The energy corrections due to finite-size 
effects arising from such a simulation are expected to be 
extremely small compared to the bulk energies considered 
here. 

The effect of TNI is important in the equation of state 
of neutron matter beyond half nuclear matter saturation 
density (p = 0.08 fm -3 ), as is clear in Fig. [TT] The two 
different TNIs added to AV8' have opposite effects: UIX 
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is repulsive, while IL7 is attractive. This is in agreement 
with the trend in neutron drops shown in Fig. 2] Our 
earlier discussion of the effects of the different terms in 
UIX and IL7 apply equally to the differences observed 
here. Furthermore, in moderately large neutron drops 
(N > 12) we have seen that the trend with JISP16 is 
similar to the trend with AV8' without TNI. We therefore 
expect that the equation of state with JISP16 will be 
similar to that of AV8' without TNIs. 

The equation of state of pure neutron matter with the 
AV8' NN interaction alone is rather soft at high densi- 
ties. It is only marginally compatible with the recently- 
observed two-solar mass neutron star [30|, HH, EH, EH- 
The relevant three-neutron force must be, in aggregate, 
repulsive at high densities [3(|. While the three-pion 
terms in IL7 (and in xPT forces) are certainly present, 
they are far too attractive in the IL7 model alone. It is 
possible to adjust the short- and long-range three neu- 
tron forces to be repulsive by varying the contributions 
of these magnitudes, as obtained in the recent studies 
of three-neutron interactions and the neutron star mass- 
radius relations [3l"| . 

VII. CONCLUSIONS 

We have computed the properties of neutron drops 
confined by external harmonic oscillator (HO) and 
Woods Saxon (WS) traps using a variety of realistic 
nucleon-nucleon and nucleon-nuclcon plus three-nucleon 
interactions (TNIs). The combination of results with HO 
and WS wells should prove useful in separating bulk and 
gradient (surface) effects, and testing the general form of 
the density functional. The pairing and spin orbit split- 
tings may also have very different behavior. 

We have employed currently available state-of-the art 
many-body methods to obtain these results and we have 
quantified the uncertainties in the results. We observe 
characteristic features such as pairing and subshell clo- 
sures in qualitative agreement with expectations. Differ- 
ences in results for the same systems are attributable in 
large part to differences in the interactions. We present 
and interpret significant sensitivity of some obervables 
to the TNI. The radii of the neutron drops appear to be 
rather robust - that is approximately independent among 



the interactions we employ. 

The results we obtain for the neutron equation of state 
as a function of density follow trends seen in the neutron 
drop results as a function of the external HO trap. This 
is significant since the neutron drops have quantified un- 
certainties on their total and internal energies as well as 
their rms radii, while it is more difficult to quantify the 
uncertainty in the calculated neutron matter equation of 
state. 

We anticipate that results for these extreme and ideal- 
ized systems may serve as guides to experiments on very 
neutron-rich nuclei. We also hope these results will in- 
form developments of improved energy-density function- 
al @,m. 
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